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Abstract. The Universe is inhomogeneous, and yet it seems to be incredibly well- 
characterised by a homogeneous relativistic model. One of the current challenges is to 
accurately characterise the properties of such a model. In this paper we explore how 
inhomogeneities may affect the overall optical properties of the Universe by quantifying 
how they can bias the redshift-distance relation in a number of toy models that mimic 
the real Universe. The models that we explore are statistically homogeneous on large 
scales. We find that the effect of inhomogeneities is of order of a few percent, which 
can be quite important in precise estimation of cosmological parameters. We discuss 
what lessons can be learned to help us tackle a more realistic inhomogeneous universe. 
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1. Introduction 

Large scale redshift surveys of galaxies reveal a hierarchy of structure that extends out 
to hundreds of Megaparsecs, yet a homogeneous cosmological model lies at the core of 
our understanding of how the Universe evolves. Such an assumption works incredibly 
well, as can be seen from current analyses of the Cosmic Microwave Background (CMB). 
The linearly perturbed Friedmann-Lemaitre-Robertson- Walker (FLRW) model can fit 
observations of the CMB with exquisite precision 

Nevertheless we cannot avoid the fact that the Universe is inhomogeneous, 
especially on small scales. Indeed, as a first approximation, space looks empty, with 
nuggets of mass and energy clustered together to form an intricate cosmic web. It is 
quite possible that the smooth description of space emerges as the average behaviour 
of a fundamentally inhomogeneous universe (cf. |2l [3l IU [5]). This is an operation that 
we are familiar with closer to home, when we describe a gas or liquid through their 
macroscopic properties instead of through a detailed atomistic description. We would 
like to be able to do the same for the Universe, i.e. to infer the smooth, large scale 
(macroscopic) properties of space time from cosmological observations without having 
to worry about the discreteness of its constituents on small scales. We attempt to 
reconcile such an observation with a homogeneous universe by saying that our local, 
granular point of view is merely the small scale limit of the description given by the 
A-Cold Dark Matter (ACDM) cosmology. We then assume that on large scales, ACDM 
is accurately described by a smooth cosmological fluid which on small scales fragments 
into particle dark matter and forms halos. Some have expressed strong concerns about 
the validity of such an approach (cf. O [71 E]), especially because it is not clear whether 
the inhomogeneities on small scales will radically affect the large scale dynamics of the 
Universe (see [9] - [27], but see also |28[ [29| |30] for counterexamples). 

There are a number of interesting and important issues that must be addressed when 
tackling this problem. In this paper we wish to focus on one: how inhomogeneities can 
affect the inferred optical properties of the Universe. This problem has been looked at 
before, throughout the 1960s [31] - [37] and 1970s [38j - [40j, and more recently [H] - 
[75] , with a variety of different assumptions and results. 

In this paper we look at how different effects can modify the redshift-distance 
relation, (z. Da)- One might expect, a priori, that studying the effect of inhomogeneities 
on {z, Da) is pointless: in 1976, Weinberg [ID] claimed that the overall {z, Da) for an 
inhomogeneous Universe should, on the whole, match that of the average homogeneous 
Universe. Views are somewhat divided on the validity of Weinberg's argument and its 
generality, with cosmologists tending to accept it wholesale, while some relativists are 
slightly more sceptical and question its assumptions (cf. [76j). In this paper we will be 
completely agnostic. 

Our approach is based entirely on the Sachs equations. Within the framework 
of these equations, we study different aspects of the redshift-distance relation. This 
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is a complementary procedure to what is usually don^. Our approach will be to 
use a number of different approximate models: high resolution N-body simulations, 
halo models, and a variety of swiss cheese models with Lemaitre-Tolman-Bondi 
inhomogeneit ies . 

While what we find is not, as yet, definitive, it does shed light on how important 
the effects of inhomogeneities might be. The paper is structured as follows. In Section 
[2] we present the Sachs equations and their various limits. In Section [3] we look at the 
effect of Ricci focusing, in Section H] we study the impact of shear, and in Section [5] we 
attempt to include local fluctuations in the expansion rate. In each of these sections, 
we need to consider different models for inhomogeneities. In Section [6] we discuss our 
results. 



2. Distance redshift relation 



The Sachs optical equations describe the evolution of optical quantities, the expansion 
of the null bundle 6 and its shear a |31j 

f^+e'+\a\' = ~R^,k-k^, (1) 

^ + 2ea = C,^^,e*"fc^e*'^fc^ (2) 

where R^js and Cap^i, are respectively the Ricci and Weyl tensors, k^^ is the null vector 
and is perpendicular to k^ and is confined to a surface tangent to a wave front. In 
the gravitational lensing nomenclature, 'shear' normally refers to the image distortion 
7. These two quantities are related but are not the same, and in order to avoid any 
confusion, we will always refer to a as to the shearing (for a detailed set of equations 
that relates these two quantities, see ||69)). 

The rate of change of the distance depends on the expansion rate of the light bundle 
in such a way that the angular diameter distance. Da-, satisfies 

^\nDA = e. (3) 

as 

Recall that the luminosity distance is given by = (1 + z^Da [13 [78]. Rewriting 
equation ([1]), we find that 

We need a set of initial conditions to solve for {z, Da), but more importantly we also 
need: (1) the Ricci curvature Ra^', (2) the Weyl curvature Cap^v] (3) the transformation 
from the affine parameter s to the redshift z. 

X The most commonly-implemented procedure to account for the effect of inhomogeneities on light 
propagation is the magnification matrix approach. In almost all cases, one further assumes that light 
propagates on unperturbed geodesies, and thus the redshift (to a source located at some fixed comoving 
distance) is calculated as in the background homogeneous models, while inhomogeneities only affect the 
niagnification/demagnification (for theoretical analysis see [42l [52l [72] , or for practical implementation 
in ray tracing codes see [49 l [50 l 173]). 



-R^pk^k^)DA. (4) 
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Using the Einstein equations {Rap — Rgap/'^ = i^T^p + ^gais) (where k = 8ttG and 
G is the gravitational constant), with the energy momentum tensor of a perfect fluid, 
one has that 

Rapk^'k^ = K(p + p)(uak''f. (5) 

In comoving coordinates, using the definition of redshift 1 + z = {uak°')e/ {Uak°')o (where 
the subscripts e and o refer to the instants of emission and observation, respectively), 
one obtains 

Rapk^'k^ = n{p + p){l + zf, (6) 

where we have used the freedom of the affine re-parametrizatioE[§| to set {uak°')o = 1. For 
the case of dust (pressureless matter), the Ricci focusing depends only on matter density 
along the past null cone. From now on, we focus only on this case (dark matter can be 
treated as being pressureless, and for the cosmological constant we have Pa + Pa = 0), 
i.e. we assume that the Ricci focusing is given by 

Rapk^'k^ = p(l + zf. (7) 

Unfortunately, there is no 'easy' trick for estimating the Weyl focusing, and so to 
calculate this factor one has to solve for the null geodesies. In fact, solving null geodesies 
is also important for linking s with z, which can be done by using the redshift formula 
together with the null geodesic equations. In comoving coordinates, we have 

l + z = {uak'')J{uak'')o = kl = ^, (8) 
which allows us to link s with z. In a homogeneous space time we have 

a = 0, 



'1 



Rapk^k^ = p^{l + z)\ (9) 

where H is the Hubble rate, and one recovers the textbook redshift-distance relation. 
If one assumes that Rajsk'^k^ = apo(l + z)^, where a = const, the Dyer-Roeder formula 
is recovered [381 1391. 



3. The role of Ricci focusing 

We will first try to quantify the effect of Ricci focusing. To do so, we need to solve (jlj) in 
terms of p{z), a{z) and s{z). In this section, we assume that a{z) and s{z) are exactly 
the same as in the background FLRW model [see the first two terms in Q]. We wish 
to study 'realistic' density profiles along the line of sight, and to do so we use density 
fields from the Millennium simulation [801 IST] . 

The Millennium simulation is an N-body simulation of the concordance cosmology. 
It consists of 10,077,696,000 particles of mass 8.6 x IO^MqH'^ within a cube of volume 

§ The afBiic parametrization is conserved with respect to hnear transformations, s — > as + 6 [79) . 
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(500/i~^ Mpc)'^. The simulation was performed using the GADGET-2 code jH2]- In our 
calculations, we use the Millennium MField, which is the dark matter density field put 
on a 256^ grid, smoothed with Gaussian kernels [Jj] of size: 1.25 Mpc, 2.5 Mpc, 5 Mpc, 
and 10 Mpc, labelled (71.25, 92.5, 95, and Qiq respectively. We randomly place an observer 
in the Millennium box, and then calculate 10^ different lines of sight (each in a different 
direction from the observer, and for each line of sight we randomly select a different 
observer). We assume periodic boundary conditions, so that when the light ray exits 
the Millennium box, it enters the other side of the box with entry angles the same as 
the exit angles (to enforce periodic boundary conditions). Since the MField consists 
of density maps on discreet time slices, to get p at any required instant we interpolate 
between different time slices. 

We calculate the light propagation out to 2; = 1.6, solve Eq. (jlj) and then write the 
distance as a deviation from the expected value in the background model (in our case 
the ACDM model), 

DAiz) = DAil + A). (10) 

Results in terms of the probability distribution function (PDF) of A are presented in 
Fig. [H We find that the larger the smoothing radius, the smaller the variance: for 
gio the standard deviation is 5.18 x 10~^ (uppermost solid curve in Fig. [1]); for g^ the 
standard deviation is 7.28 x 10~^ (the solid curve second from the top in Fig. [1]); for (72.5 
the standard deviation is 9.52 x 10^^ (the solid curve third from the top in Fig. [T]); and 
for (71.25 the standard deviation is 1.21 x 10^^. It is interesting to compare these results 
with the ones obtained using the weak lensing formula (in the Born approximation), 

Xe 

Awl = —H'^Qm [ dx^^xa-'S{x), (11) 

J Xe 


where x is the comoving coordinate, dx = dz/H{z) and we use the same S{z) as before. 
The results are plotted with dotted lines in Fig. [H As shown, the PDFs are very similar 
to the ones obtained within the Ricci focusing regime. Finally, we calculate the standard 
deviation by squaring the expression above (where the mean is zero by construction) 
and replacing 6"^ with the matter power spectrum [83], 



In this case, the standard deviation of the distance correction i^ 1.45 x 10~^ (the 
dashed curve in Fig. [1] is the Gaussian PDF with mean of zero and standard deviation 
of 1.45 X 10-2). 

II P = ^iW, where W ~ ex.p[x^ /a'^]. 

% Here, in order to compare with the Millennium simulation, we use the cosmological parameters that 
were used in the Millennium simulation. If the WMAP7 set of cosmological parameters is used instead, 
then CTA = 1.57 x 10"^ 



'1 I \Xe X 

.1 + 2; X 

Xe 



dk k 



P{k,z) 
2n 



(12) 
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Figure 1. PDF of A. Solid curves (from tlie top) are tire results obtained by solving 
Eq. 13] within the Ricci focusing regime, for the Millennium maps 1710, (75, 52.51 
51.25- The dotted lines are the results obtained from (fTTj) for the same density fields. 
The dashed line is the Gaussian PDF with the variance calculated from Eq. (fT2|l . 



The above results clearly reveal a pattern. First of all, the smaller the smoothing 
radius, the larger the variance, although even with the 1.25 Mpc smoothing radius 
we get a smaller standard deviation than was obtained within the framework of linear 
approximation. Also, note that the maximum of the PDF lies on the demagnification 
side (i.e. A > 0), so a random object is most likely to be dimmer than on average. 
Nevertheless, in all cases the mean is almost zero: 1.20 x 10^^, 1.38 x 10"^, 1.29 x 10~^, 
1.44 X 10"^, for 5(1.25, 92.5, 95, and 910 respectiveljH 



4. The role of smoothing and Weyl focusing 

We now extend our analysis to investigate the role of smoothing and the effect of Weyl 
focusing on the optical properties of the Universe. When dealing with the Weyl focusing, 
it is not straightforward to link density fields with the Weyl curvature. An ideal situation 
would be to have an exact generic solution of the Einstein field equations, and use that to 

+ Note that, for each different line of sight, we chose a different observer location. If we located 
an observer inside some deep inhomogeneity instead, such as a void of ^1.25 = —0.84, ^2.5 — —0.78, 
^5 -0.66, 610 = -0.43, then the mean would be 2.54 x 10"*, 2.58 x lO"'*, 2.46 x lO""*, 2.07 x lO""^. 
This suggest that the effect of local inhomogeneity on the distance (to z = 1.6 in our case) is of order 
10"''. In this paper we try to avoid the 'local' effects introduced by nearby inhomogeneities, so that we 
have a more clear insight to light propagation effects alone. For papers that study the local bias see 
[84]-[89]. 
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Figure 2. The thick sohd Hne is the PDF obtained in the full halo model. The dashed 
line is a variation of the halo model, where we switched off the Weyl focusing. The 
dotted lines (from the top) are halo models with A so equal to 80 and 40 respectively 
(see text for explanation). For comparison, we also present the results when the gi.25 
Millennium map is used (dot-dashed line). 



model the evolution and effects of the Weyl curvature. Since there is no such a solution, 
in practise we must always employ some approximations. Here we will evaluate the 
Weyl focusing within the halo model. 

Simply put, the halo model considered here consists of a universe filled with different 
halos of masses above 10^^ H'^Mq, while the rest of the mass density (in objects of masses 
below this threshold) is assumed to be distributed homogeneously. The halo model we 
use closely follows the construction outlined in detail in Ref. [64j . Each halo is described 
by the Navarro- Frenk- White (NFW) density profile |90] . 

P^'^=P-^'\r/R)illr/Rr ^''^ 
where R is the radius of the halo, defined as 



M 



1/3 



^ [{4/3)7rpUz)Aso) ' ^^^^ 
with Pm{z) = po{l + zY, and Aso = 180 (the halo is defined as an overdensity of mass 
M, whose constant density contrast is 180 with respect to the mean matter density, cf. 
[9TI l6l] ). Integrating p(r), and using the expression above, one finds that 



6r 



A 



so 



3 ln(l + c) -c/(l + c)' 



(15) 
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where c is the concentration parameter, whose evolution given in [22] as 

c = 10.14[M/(2 X lO^2/,-i;vf0)]-°-°si(l + (16) 
The number density of halos at a given instant in time is given by the halo mass function 

with / and a given by [91] 

f{a,z) = 0.301exp[-|lna"i + 0.67|^-^2|], 

oo 

a{M,z) = / k'^P{k)W^{k,M)dk, (17) 

271^ J 


where G{z) is the growth factor (this is usually written as D, but we use Q here in order 
to avoid confusion with distance). 

The halo model allows us to address the problem of smoothing — instead of dealing 
with a continuous density field, we now have a discrete set of dark matter halos — 
although one should keep in mind that, given the characteristics described above, this 
is not really a completely 'discrete' model. It also allows us to take into account the 
effect of the shearing, by assuming that the Weyl focusing for a particular halo, i, is 
calculated using the Lemaitre-Tolman (LT) model [93l [94] (an exact general relativistic 
model of an spherically symmetric, inhomogeneous, non-stationary space time) [171 HE] 

^ C^p,.e*'^k^e*^k^ = 1^ (p _ p) , (18) 
where b is the impact parameter, R is the areal radius. 

We apply the weak field approximations (cf. Ref. [33j) and assume that the total 
Weyl focusing is a sum over all contributors, C = Note that, outside the halo, 

Ci ~ b'^/R^, and so only halos that are close to the light ray contribute most significantly. 

We again solve Eq. for 10^ different lines of sight. The Ricci and Weyl focusing 
are calculated using the halo model, while s{z) is assumed to be the same as in the 
FLRW background model [see the first term of ([9])]. The resulting distribution for A is 
presented in Fig. [2] (thick solid curve). As expected, the variance is much larger than 
when the smoothed Millennium maps were used. Also, the maximum is shifted more 
towards positive values of A. The mean, as before, is still negligible, at 8 x 10~^. 

For comparison, and to gain a better understanding of the different factors, we now 
consider rather unphysical variations of the halo model. The first variation is when we 
neglect the shearing - when solving we put a = 0. The result of this is shown as 
the dashed line in Fig. [2l As can be seen, it is significantly shifted to the magnification 
side (positive A) - the mean is 1.02 x 10^^. This shows that, when considering the halo 
model, the role of the shearing should not be neglected when calculating the distance 
correction. 
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To better understand the relation between the halo model and the smoothed 
Millennium density fields, we can study two more unconventional modifications of the 
halo model. In order to reduce the level of discreteness of the halo model, we decrease 
the value of the parameter Aso, which results in larger halos and a lower amplitude of 
the halo density profiles. In the first case we reduce Aso to 80, which means that the 
radius of the halo increases by roughly 30%. In the second case, we reduce Aso to 40, 
which results in a 65% increase in the halo's radiufl We plot these results as dotted 
lines in Fig. O As can be seen, when we increase the halo radii (and hence decrease 
the degree of discreteness), we approach the results obtained when using the smoothed 
Millennium maps. 

So far, the results reveal the same pattern: a skewed distribution with a maximum 
at the demagnification side {Amax > 0). This can easily be understood within the 
framework of the weak lensing formula ( ITTl) : when light rays pass through the large- 
scale structure, they are more likely to propagate through voids than through more 
dense, compact structures. Therefore, it is more likely to have 6 < along the line of 
sight and, as follows from (fTT]) . 5 < ^ A > 0. This is clearly visible when we use the 
smoothed Millennium maps. If the smoothing scale is sufficiently large, the resulting 
density field approaches a Gaussian field (all non-linearities are washed away, and we 
recover symmetry between overdense and underdense regions), and so the maximum 
approaches zero (for pure Gaussian fiuctuations Amax = 0). In the halo model, we deal 
with halos that are very compact and occupy a relatively small volume, which makes 
propagation through voids even more likely to occur. We therefore see a shift of Amax 
towards higher demagnification (see Fig. [2]). Occasionally, light rays pass close to a 
halo, and then both the Ricci and Weyl focusing magnify the light bundle. Thus, on 
average, the total balance is recovered, and the mean of A is almost zero (in the lensing 
framework, as seen in f lTT]) . due to matter conservation it is exactly zero). Here we find 
a mean of order 10~^. A small deviation from the expected value in the FLRW model 
is expected. For example the analogue of the Integrated Sachs- Wolfe effect will affect 
the redshift, and thus {z, Da)- Also, the local environment plays a role - note that our 
observer is located inside an inhomogeneity (i.e. not within a homogeneous region) - 
see footnote on page 6. 

The same pattern has also been reported in previous studies - see for example in 
|42] . and in ray tracing though N-body simulations. In fact, our halo model provides very 
good agreement with the ray tracing through the Millennium simulation, cf. ^91 150] . 
It is also in almost perfect agreement with the PDF generated using the turboGL 
code|59l IM jjl . 

* For example, a halo of mass lO^^A/© and A50 = 180 has R « 330 kpc [this follows from ([T4| ]. for 
^so — 80 and A50 — 40, R « 430 kpc and R ~ 545 kpc respectively. 

tt The turboGL code uses the halo model to calculate the PDF of the weak lensing convergence. The 
effect of the shear is neglected in turboGL. However, if the zero-shear PDFs are compared, they are 
in agreement. Note that, as seen in Fig. [21 the zero-shear PDF, apart from a small change in its 
amplitude, looks like the translation of the nonzero-shear PDF. 
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The results are therefore consistent with the most commonly-implemented 
procedures, where the effects of inhomogeneities on light propagation are calculated 
using the magnification matrix approach. In almost all cases, one further assumes that 
light propagates on unperturbed geodesies, and so the redshift (to a source located at 
some fixed comoving distance) is calculated as in the background homogeneous models, 
while inhomogeneities only affect the magnification/demagnification. So far, we have 
been doing the same (that is, we assumed that s{z) is the same as in the background 
FLRW model). We therefore examine the effect of changing s{z) in the next section. 

5. The role of the non-uniform expansion rate 

The last effect that remains to be examined is related to s{z). So far we have been 
applying the FLRW formula, but as shown in ([8]), to calculate it correctly we need to 
solve null geodesies. We now address this point by studying light propagation within 
the LT Swiss Cheese model. The advantage of the LT model is that it is an exact model, 
and so no approximations are required when modelling the light propagation, evolution 
of perturbations, or the Weyl curvature. The disadvantage is that each particular 
inhomogeneity must be spherically symmetric, and we do not have the freedom to mimic 
the Millennium density profile. Furthermore, we cannot describe virialised objects such 
as massive matter halos (we will come back to this point later on). 

Let us first consider a 'mild' Swiss Cheese model, i.e. one without large density 
fluctuations. Each inhomogeneity has a radius of 20 Mpc, beyond which the system 
becomes homogeneous. The density contrast inside the inhomogeneity is generated 
using the log-normal PDI^, with cr^ = 0.96 (i.e. R = 10 Mpc). An example of a 
random line of sight through this model is presented in the panel to the lower left of 
Fig. 131 Within this model, we solve the null geodesic equations using ([8]) to obtain s{z). 
The Weyl focusing is calculated from f lTSj) and the distance is calculated by solving 
In order to avoid effects of local structures, we place the observer and the source within 
the homogeneous regions. The resulting PDF of A for a source at z = 1.6 is presented in 
Fig. m (dotted line). The mean is 5.75 x 10~^. Qualitatively, the results are very similar 
to the ones obtained before: the PDF is skewed, with the maximum at the positive side 
of A. Quantitatively, however, the variance is much smaller than in previous cases. This 
is because of the smaller amplitude of inhomogeneities along the line of sight, which can 
be seen by comparing the upper left (Millennium) and upper right (halo model) panels 
of Fig. [3] to the lower left panel of Fig. |3] (the 'mild' Swiss Cheese case). 

To bring us closer to the Millennium case, but in such a way as to include shearing 
and a modified s{z) relation, we now consider a more 'realistic' model, with large deep 
voids (with, for example, a present day radius of 14 Mpc and a density contrast of —0.83). 
These voids have highly non-linear walls which, due to the high density contrast, are 

jiy, where cr^; = ln[l + aj^}, and aj^ = 

^ J dkV{k)W^{kR)k'^. 





t PiS) 



: exp 



(ln(l+a)+o-^,/2)^ 
2^, 
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Figure 3. Upper left: Density along a random line of sight through the 51.25 density 
field of the Millennium simulation. Upper right: Density along a random line of sight 
through the halo model. Lower left: 'Mild' Swiss Cheese model. Lower right: Density 
along a random line of sight through the highly non-linear Swiss Cheese model (solid 
line - 'u'-shaped pattern) and the contrast of the expansion rate Sh — H/H—1 (dotted 
line - note that in regions where S is of high amplitude, Sh is negative). 



at the collapsing stage, clearly visible in the lower right panel of Fig. |3] (dotted line). 
The PDF obtained within this model is presented in Fig. H] (solid line). The mean is 
— 1.01 X 10~^. The first qualitative difference is that the maximum is now at negative 
A (magnification). Secondly, due to large fluctuations in the density field, the variance 
is larger, and is comparable with the giQ Millennium density field (dotted-dashed line 
in Fig. H]). 

The surprising result is that the maximum is on the magnification side {Amax < 0). 
Such a phenomenon has not been reported before. Let us closely examine what could 
lead to such an unexpected behaviour. Three working hypotheses can be put forward to 
explain this: (1) the effect of the Weyl focusing, (2) the problem with collapsing regions, 
(3) an artefact of the symmetry of the density contrast. With regards to the shearing, 
we have seen in the previous section, when comparing the halo model with and without 
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Figure 4. PDF of A in various models. The central dotted line is the PDF obtained 
within the 'mild' Swiss Cheese model. The PDF of the highly non-linear LT Swiss 
Cheese model is given by the solid line. The dashed line is also the PDF obtained in 
the non-linear LT Swiss Cheese model, but when we neglect the effect of collapse on 
the redshift relation (see text for details). For comparison, the PDF for the 510 density 
map (from Fig. [1]) is shown as the dashed-dotted line. 

the Weyl focusing, that if cr 7^ 0, the PDF shifts to lower values of A. However, as 
seen from Fig. |5l the role of the Weyl focusing for our Swiss Cheese models is small 
(cf. [To] and [72] for a similar conclusion regarding the role of the Weyl focusing in this 
type of model). Essentially, if the density contrast is not high (unlike in the halo model; 
compare the upper right panel of Fig. [3] with the lower right panel), the role of the Weyl 
focusing is small (compare the set of 3 solid lines with the set of 3 dashed or dotted 
lines in Fig. [5]). We therefore conclude that the shear cannot be responsible for the 
shift of Amax from demagnification to the magnification side. This can be confirmed by 
re-running the same model with a set to zero - the resulting PDF is exactly the same 
as in the case of a 7^ (solid line in Fig. H]). 

Our second hypothesis involves considering the role of collapsing regions. In this 
case, H < 0, which may strongly affect s{z). In the real Universe, whenever the 
density field has a high amplitude, cosmic structures are virialised, and hence they 
are not collapsing. Unfortunately, there is no rotation within the LT model, which 
could prevent the collapse. Therefore, as seen in the lower right panel of Fig. |3l the 
expansion is negative within the wall. In order to 'correct' the model for the virialisation 
effect, whenever the expansion rate is negative, we assume that the structure is already 
virialised and so the redshift does not change within this region. That is, we assume 
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that d/c^/ds = in these regions, whereas in all other regions we use the exact formula, 
i.e. dfco/ds = -RR{k'^)^ - R'{R'/{1 + 2E)){k^f . We should exercise caution when 
analysing these results, as this 'correction' was done by hand. The fact that some 
regions are virialised, and hence d/c^/ds = 0, is not a consequence of a model, and so 
there is clearly an issue of self-consistency. The resulting PDF is presented using a 
dashed line in Fig. |U where we find that the maximum of the curve is shifted even 
further towards negative values of A. The mean is almost twice as large as before, 
—1.91 X 10^'^. Clearly this does not explain why A^ax < 0, and in fact it makes things 
even 'worse'. 

Our final hypothesis concerns the symmetry of the density perturbations. We can 
already see that this cannot be the cause, as A^ax > for the 'mild' Swiss Cheese 
model. Nevertheless, we return to the Millennium smoothed maps and perform further 
tests. We proceed as before, with the only difference being that, when solving s{z), we 
no longer assume a uniform expansion rate as in (Q, but instead assume that 

^^ = {l + zyH{l + 6H), (20) 
where H is the expansion rate of the background model. In general 

HSh = ^ (e - 6) + aabe'^e' 

where G is the matter scalar of the expansion and aab is the shear of the matter velocity 
field. We further work under the assumption of zero shear and so we only consider 
fluctuations in the expansion field, which in the linear regime are given by 

Sh = —fS. (21) 

Here, / is the growth function, / = Q/{QH). This method has thus 2 limitations, and 
works well as long as (1) the density contrast is not too high and (2) the matter shear 
is negligible. We proceed further with this toy model, as it will bring more insight into 
results obtained within the exact LT model. 

Naively, one can think of this as follows: because (5) = = {5h)i we should expect 
(qualitatively) the same results as before, except for a potentially larger variance due 
to the presence of both 5 and 5h- Our results are presented in Fig. [6l First of all, the 
variance changes only slightly. Secondly, and most importantly, we find that A^az < 0. 
In the giQ case (with 10 Mpc smoothing scale), the results are strongly comparable with 
our non-linear Swiss Cheese model; for comparison we present its PDF in the upper left 
panel of Fig. [6] (dot-dashed line). Also, the mean of A is of the same order of magnitude 
as in the Swiss Cheese model, i.e. ~ 10~^ — see Table [H 

Finally, as before, we try to 'correct' the model for the effects of virialization. As 
seen from f l2T]) . if 

. . 3 
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Table 1. The mean and standard deviation of A for a source at z — 1.6. 



MnHpl 


— 

UH — u 


UH 7= u 
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-8.77 X 10- 
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-2.85 X 10" 
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-3 
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non-linear Swiss Cheese 




-1.01 X 10- 


-3 


-1.91 X 10- 


-3 


halo 


-8.02 X 10"^ 






-2.31 X 10" 


-2 



Table 2. Middle column: The percentage of regions within the Millennium simulation 
whose density contrast is higher than the threshold above which J^f < — 1 Right column: 
The average expansion rate after setting = — 1 wherever Sh < —1- 

Map % of 5 > 6t {6h)t 



c/1.25 1.459% 1.95 X 10-2 

^2.5 0.762% 5.31 X 10-3 

g5 0.089% 2.69 x 10"^ 

910 



then such a region collapses. In the (71.25 case, almost 1.5% of regions have a density 
contrast above this threshold (see Table [2]). As before, we set 

if 5 > 5t =^ = -1- 

In this case, if we calculate the average 6h, due to the fact that we neglect the most 
negative contributions, the average is no longer zero, and in the case of (71.25 the average 
5h is approximately 0.0211 (it can only be zero if 6h < —1)- The resulting PDFs are 
presented in Fig. [6] (dashed lines). As in the case of the non-linear Swiss Cheese model 
before, the PDF shifts towards even more negative values of A. For detailed values of 
the means of the PDFs, see Table [H 

Proceeding further with the idea of a variable expansion rate and its impact on 
s{z), we also recalculate our results for the halo models that were considered previously. 
Since halos are virialised, their interiors are set to 6h = —1- The resulting PDF is 
presented in the right panel of Fig. [71 As was the case for the Millennium maps, a shift 
towards magnification is evident. This should be compared with the mild Swiss Cheese 
model, or with the models which assumed a uniform expansion rate, where the PDF 
peaked on the positive side of A. 

In Sees. |3]and|l]we assumed that s{z) was the same as in the FLRW background 
model. As a result, we found that (A) ?a 10"'^. For the Swiss Cheese model, we had 
variations in s{z) with regards to the background, but still found that (A) ^ 10"^. 
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Figure 5. The Weyl focusing (i.e. the ratio of the shearing to present day matter 
distribution cr^/po) in the different models that were considered. From bottom: the 
three dotted hnes are 3 examples of the shearing in the mild LT Swiss Cheese model. 
The shearing at 2; = 1 is of order 10~^ — 10"''. Above those, the next 3 (dashed) 
lines are 3 examples of shearing in the highly non-linear Swiss Cheese model. Here, 
at z = 1, the shearing is of order 10"^. The next 3 (solid) lines show shearing in 3 
examples of the halo model. At z = 1, the shearing is of order 10~^. For comparison, 
the uppermost dash-dotted line is (1 -I- z)^ (the evolution of matter density). 



It was only when we started to 'correct' our models for virialization that we found 
deviations of order of 10~^ (for the halo and (yfi.25 models). We now want to follow 
this up further, and test how the average of the distance follows the background's s{z). 
Note that when we 'force' some regions not to collapse (by setting 5// = — 1 whenever 
5// < — 1), we change the average expansion history (see Table |2] and Fig. [8]), and thus 
also s{z). This leads us to question: if we recover the average s{z) correctly, will we 
obtain the correct mean of the distance, without the need for a detailed calculation of 
the Ricci and Weyl focusing? 

We study this issue within the halo and (71.25 models, as they have the largest 
deviation of (A). The PDFs obtained within these models are presented in Fig. [71 
The mean of A is shown as a vertical solid line. The average of the expansion rate 
fluctuations calculated on a slice of constant z is presented in Fig. [HI Before proceeding 
further, we must comment on one important aspect of our computations. Cosmologists 
tend to use the redshift as a proxy for time, so they express the evolution of different 
fields in terms of the redshift, for example: 



p{s)^ p{z) = Poil + zf, 
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Figure 6. Upper left: Millennium map giQ with 5h = ^ (dotted line, cf. Fig. [T]), 
(5^ ^ (solid line), and (5// > — 1 (dashed line, in this case indistinguishable from the 
solid line). For comparison, the PDF of the non-linear Swiss Cheese model is shown 
as the dash-dotted line. Upper right: Millennium map 35 with 5h = ^ (dotted line), 
5fj ^ Q (solid line), and (5^f > — 1 (dashed line, hardly distinguishable from the solid 
line). Lower left: Millennium map (72.5 with 5h — ^ (dotted line), 5h (solid line), 
and Sh > —1 (dashed line). Lower right: Millennium map 51.25 with Sh ~ (dotted 
line), 6h ^ (solid line), and 6h > (dashed line). 



H{s) ^ H{z) = Ho^/amii + + + zy + n^. (22) 

Only in the case of perfect homogeneity do we have that z (the proxy for time, i.e. the 
redshift of a photon at a given time in a perfect homogeneous model) coincides with z 
(the actual shift of the photon's frequency). Thus, in solving the Sachs equations, one 
must solve the following: 

dz 

- = (1 + zf{H) = (1 + zfH{z){l + 5h{z)) 
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Figure 7. PDF of A. Solid line: mean. Dotted line: solution of (|23p . Dashed line: 
when (j22p is used together with the assumption that z = z. Left: halo model, Right: 
91.25 Millennium map. 
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Figure 8. Average of the fluctuation of the expansion rate at constant z. Solid line 
51.25 Millennium map, dashed line for the halo model. 



^ = (1 + 5^5). (23) 

We solve these, and plot the results as dotted lines in Fig. [71 As can be seen, the 
difference between the actual value, and that obtained from the above set of equations, 
is of order of 10~^, (as before, where we set s{z) to its value in the background ACDM 
model). For comparison, we solve the above equation with 

z = z R^pk'^k^ = po(l + zY- 

This time, the result (given by the dashed line) is significantly displaced to the left. The 
reason for the shift towards higher magnification is that if {6h) > then, as follows from 
(123|) . we have that z > z, and so if we replace z with z, we increase the Ricci focusing, 
hence higher magnification. This simple exercise shows that if one is interested in having 
an average Da{z) without having to worry about the Ricci or Weyl focusing, one can 
simply use the average s{z) of one's model. Caution must be taken, though, when 
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studying models in which p is not proportional to {1 + z)^. Such models include, for 
example, Gpc-scale inhomogeneous LT model^ and some backreaction models. 

In this section we have tested all 3 factors that affect the relation between z and Da 
via the Sachs equation: p{z), cr{z), and s{z). In order to model them in a self-consistent 
manner, we used the LT Swiss Cheese model. In all cases, we found that the distance 
does not change, on average, by more than 10"'^ (at z = 1.6). This result is fairly 
consistent with what has been found by other groups who have employed exact models, 
and ensured the proper randomization of light rayj§| (cf. [171 H8l l56| \7Q\ 1711 175]). 

Our main motivation has been to follow, as closely as possible, the standard 
procedures to test self-consistency and the impact of inhomogeneities on the optical 
properties of the Universe within models that follow the FLRW evolution. Nevertheless, 
it is not completely obvious that one should recover the FLRW results for these models, 
unless one completely agrees with Weinberg's reasoning. If one does take Weinberg's 
results at face value, then our 10^^ results should be considered as a big deviation, about 
two orders of magnitude larger than what is allowed by the analogue of the Integrated 
Sachs Wolfe effect (which should be of order 10~^). 

From one point of view, our results have turned out to be quite surprising. In the 
case of large fluctuations in the density and expansion fields, the mode of A shifts to the 
magnification side. Such a behaviour has not been reported before. The effect is only 
visible if the degree of inhomogeneity is large - for small amplitude Swiss Cheese models, 
we recover the standard results. This can be explained using the same 'voids argument' 
that we used in the case of the weak lensing approximation - since the structures in the 
Universe form a cosmic web, it is more likely for photons to propagate through voids 
(hence the position of the mode is associated with the amount of void regions). As seen 
from (fTTj) . the weak lensing regime is only sensitive to density fluctuations, and so the 
mode is positive. If we introduce fluctuations in the expansion rate, then the redshift 
decreases, as seen from (l20l) . Hence an object located at the same distance has a higher 
redshift, and thus gives the impression of being magnified compared to the homogeneous 
cas^. This is the reason why the mode is now on the magnification side. 

A valid point, and one worth further investigation, is related to the fact that, due 
to the virialisation of high density regions, the average expansion rate may deviate 
from the background value. This was already pointed out in [98| [99] in the context 

I In fact the results of this paper provide a fresh view into different configurations of the LT model, 
which consist of giant voids |95[|96] and giant humps ^7\. Within the giant voids p increases outwards, 
which results in the decrease of the expansion rate along the past null cone, hence the dimming of 
the supernova within the giant voids is caused by the change of the s{z) relation. Whereas for the 
giant hump models p decreases outwards making the Ricci focusing less efficient, and hence causing 
additional dimming. 

§ Note that for models where the distance correction is calculated using the Born approximation, have 
by the construction (A) = [52 l [59 l [M l [72] 

II The same fact can also be visualised in terms of the slope of the relation (z, Da) - when the redshift 
increases faster than in the homogeneous case, the slope decreases. See, for example, Fig. 4 of [74], 
which shows the change of the slope in the relation {z,Da), depending on fluctuations in the density 
and expansion rate. 
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of backreaction models, although as shown in |lUm IIUI) . the presence of matter shear 
may decrease the backreaction effect, leaving the evolution of the average quantities 
relatively unchanged. Nevertheless, in such cases, the 3D volume average on surfaces of 
constant time may not coincide with the average on slices of constant z. As a result, as 
we have shown above, we would still expect a change in the redshift-distance relation. 
However, we reiterate that care should be taken in interpreting these results - changing 
H by hand may not lead to self consistent results (although keeping the expansion rate 
uniform may also turn out to be inconsistent). To test this phenomenon more closely, 
we require exact models that allow for virialisation. The LT models are definitely not 
suitable for such a study. 

6. Conclusions 

The overarching question that we wish to address is: does the FLRW model, which 
correctly describes the evolution of the background, correctly recover (on average) the 
distance-redshift relation? Our results, while not entirely conclusive, are suggestive. 
We have explored a suite of models which are inhomogeneous, but close to the FLRW. 
We find that, if the background model correctly describes the evolution of the average 
expansion rate and the average matter density, then indeed the homogeneous model is 
a good approximation. 

Our results show that Ricci focusing alone introduces deviations from the average 
value. The PDF is skewed and the maximum is on the demagnification side (a random 
object will most likely be dimmer than the average). The change of the mean is 
negligible, which means that, on average, distances coincide with the expected distance 
in the background model. The deviation from the average (for a particular line of sight) 
depends crucially on the smoothing scale and, unsurprisingly, the smaller the smoothing 
scale, the larger the deviation. We find that Weyl focusing is important for objects with 
very high density contrast, and yet even in this case, the bias in the mean of A is 
negligibly small. 

We also find that the presence of large fluctuations in the expansion rate also 
affects the final results: there is a slight shift of the maximum of the PDF towards the 
magnification side (a random object will most likely be brighter than the average), and 
the variance changes. In particular, we find that if the average expansion rate deviates 
from the background expansion rate, the distance (on average) changes by a few percent. 
In the real Universe, such a situation may occur in the late non-linear stage of evolution, 
where large voids expand much faster than the background, while virialised over-dense 
regions do not expand at all (for these fluctuations to cancel each other, at some stage 
the presence of collapsing regions is required). 

In summary, we have studied the various biases that cosmological inhomogeneities 
can introduce into the inferred optical properties of the Universe by considering a range 
of toy models. There are, of course, limitations to the accuracy and self-consistency of 
these models, and how well they mimic the real Universe. We find that the presence of 
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inhomogeneities (around a fixed background) changes distances by at most a few percent, 
altliougli tlie average is closely related to the background model. The only possible 
deviation of the average of the distance from the distance predicted by the background 
model is in situations where the background model does not correctly describe the 
evolution of the average quantities. 

The lessons we have learned can be summarized in the form of a message to our 
readers. For cosmologists, we would like to point out that, with the increasing precision 
of cosmological observations, a change of few percent might soon become important for 
accurate estimation of cosmological parameters (cf. |63 | 167 ) ). In particular, the PDFs 
are skewed and hence, for a limited number of observational data, there might be a 
bias due to the fact that the mode does not coincide with the mean, or median. For 
relativists, we would argue that, to study the effect of inhomogeneities on observations, 
it is first important to understand how much inhomogeneities affect the evolution of 
the background model. Such a study should focus on p{z) and H{z), and especially 
on potential deviations from f l22l) . as this has the biggest impact on the distance- 
redshift relation (see, for example. Fig. [7] for a few percent change in the case where 
p{z) ^ (1+2;)^ (cf. [68]). Thus, pinning down p{z) and H{z) will also provide the average 
of Da{z), without the need to account for the Weyl focusing or to model complex aspects 
of the Ricci focusing). 
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